#ifndef _HFSOLVER_K_H_
#define _HFSOLVER_K_H_
#include<iomanip>
#include<algorithm>//for sort algorithm
#include<cmath>
#include<eigen3/Eigen/Eigen>//matrix library

#include"orbital.h"
#include"system_k.h"
#include"help.h"

using Eigen::Matrix;
using Eigen::Dynamic;
using Eigen::SelfAdjointEigenSolver;

using Eigen::ComplexEigenSolver;
using std::cout;
using std::stable_sort;
using std::abs;
using std::setw;

///spherical HF solver,for System_WS
/* class SHFSolver_WS */
/* { */
/*  public: */
/*   typedef System_Tz_WS System_TzType; */
/*   typedef typename System_TzType::DataType DataType; */
/*   typedef typename System_TzType::OrbitalType OrbitalType; */
/*   typedef Matrix<DataType,Dynamic,Dynamic> Mat; */
/*   typedef Matrix<DataType,1,Dynamic> Vec; */

/*  SHFSolver_WS(const System_TzType * Ptr):pSystem_Tz(Ptr){} */
/*  SHFSolver_WS(const System_TzType * Ptr,int _A,int _Z):pSystem_Tz(Ptr),A(_A),Z(_Z){} */

/*   void setAZ(int _A,int _Z){A=_A;Z=_Z;} */

/*   ///return onebody part of the hamiltonian */
/*   DataType getOneBody(int bra,int ket) */
/*   { */
/*     return pSystem_Tz->get1B(A,bra,ket); */
/*   } */
/*   ///return twobody part of the hamiltonian <ab->J|V|cd->J>, */
/*   /\*noting |ab->J> means sqrt(1+delta{ab})|ab,J>*\/ */
/*   DataType getTwoBody(int a,int b,int c,int d,int J) */
/*   { */
/*     return pSystem_Tz->get2B(A,a,b,c,d,J); */
/*   } */
/*   ///1/(2j+1) sum_h{ (2J+1) < bra h->J | V2B | ket h->J >}, h is hole */
/*   DataType SHFDiagram(int bra,int ket) */
/*   { */
/*     int jj=pSystem_Tz->Orbitals[bra].jj; */
/*     DataType val2B=0; */
/*     for(int i=0;i<HFOrbitals.size();i++) */
/*       { */
/* 	if(HFOrbitals[i].state) continue;//skip particles */
/* 	int Jmin=abs(jj - HFOrbitals[i].jj)/2; */
/* 	int Jmax=(jj + HFOrbitals[i].jj)/2; */


/* 	int GroupIndex=HFOrbitals[i].GroupIndex; */
/* 	int dim=pSystem_Tz->Groups[GroupIndex].size(); */

/* 	for(int k=0;k<dim;k++) */
/* 	  for(int l=0;l<dim;l++) */
/* 	    { */
/* 	      for(int J=Jmin;J<=Jmax;J++) */
/* 		{ */
/* 		  int brak=pSystem_Tz->Groups[GroupIndex][k]; */
/* 		  int ketl=pSystem_Tz->Groups[GroupIndex][l]; */
/* 		  val2B += getTwoBody(bra,brak,ket,ketl,J) * HFOrbitals[i].Coeff[k] * HFOrbitals[i].Coeff[l] * HFOrbitals[i].OccuP * DataType(2*J+1); */
/* 		} */
/* 	    } */
/*       } */
/*     val2B/=(jj+1); */
/*     return val2B; */
/*   } */

/*   ///HF single particle hamiltonian, <bra|V1B|ket> + 1/(2j+1) sum_h{ (2J+1) < bra h->J | V2B | ket h->J >}, h is hole */
/*   ///(group,a) corresponding bra,(group,b) corresponding ket */
/*   DataType HFhamil(int group,int a,int b) */
/*   { */
/*     int bra=pSystem_Tz->Groups[group][a]; */
/*     int ket=pSystem_Tz->Groups[group][b]; */
/*     return getOneBody(bra,ket)+SHFDiagram(bra,ket); */
/*   } */
  
/*   void printOrbitals(int num=-1) */
/*   { */
/*     if(num==-1) */
/*       num=HFOrbitals.size(); */
/*     for(int i=0;i<num;i++) */
/*       { */
/*       	cout<<HFOrbitals[i].l<<"\t"<<HFOrbitals[i].jj<<"\t"<<HFOrbitals[i].tz<<"\t"<<HFOrbitals[i].e<<"\t"<<HFOrbitals[i].OccuP<<endl; */
/*       } */
/*   } */
  
/*   ///initialize */
/*   void init(); */

/*   ///sort and fill orbitals. */
/*   void fillOrbitals(); */

/*   ///construct matrix and diagonalize, generate new HF orbitals */
/*   void diag(); */
  
/*   ///iteration */
/*   void iter(); */

/*   ///cal. HF energy */
/*   void calEnergy(); */
  
/*   int A,Z; */
/*   vector< SHF_Orbital<DataType> > HFOrbitals; */
/*   DataType Energy; */
/*  private: */
/*   const System_TzType * pSystem_Tz; */
/* }; */


//===============================================================================================

///spherical HF solver,for System_HO_Tz_K
class SHFSolver_K
{
 public:
  typedef System_Tz_K System_TzType;
  typedef typename System_TzType::DataType DataType;
  typedef typename System_TzType::OrbitalType OrbitalType;
  typedef Matrix<DataType,Dynamic,Dynamic> Mat;
  typedef Matrix<DataType,1,Dynamic> Vec;

 SHFSolver_K(const System_TzType * Ptr):pSystem_Tz(Ptr),numIter(0){}
 SHFSolver_K(const System_TzType * Ptr,int _A,int _Z):pSystem_Tz(Ptr),A(_A),Z(_Z),numIter(0){}

  void setAZ(int _A,int _Z){A=_A;Z=_Z;}

  ///return onebody part of the hamiltonian
  DataType getOneBody(int bra,int ket) const
  {
    return pSystem_Tz->get1B(A,bra,ket);
  }
  ///return twobody part of the hamiltonian <ab->J|V|cd->J>,
  /*noting |ab->J> means sqrt(1+delta{ab})|ab,J>*/
  DataType getTwoBody(int a,int b,int c,int d,int J)
  {
    DataType v;
    if(numIter==0)
      {
    	v=pSystem_Tz->get2B(a,b,c,d,J);
      }
    else if(numIter==1)
      {
    	v=pSystem_Tz->get2B(a,b,c,d,J);
    	TwoBodyMat.push_back(v);
      }
    else
      {
    	v=TwoBodyMat[count];
    	count++;
      }
    return v;
  }
  ///1/(2j+1) sum_h{ (2J+1) < bra h->J | V2B | ket h->J >}, h is hole
  DataType SHFDiagram(int group,int a,int c);


  ///HF single particle hamiltonian, <bra|V1B|ket> + 1/(2j+1) sum_h{ (2J+1) < bra h->J | V2B | ket h->J >}, h is hole
  ///(group,a) corresponding bra,(group,b) corresponding ket
  DataType HFhamil(int group,int a,int b)
  {
    int bra=pSystem_Tz->Groups[group][a];
    int ket=pSystem_Tz->Groups[group][b];
    double t1=getOneBody(bra,ket);
    double t2=SHFDiagram(group,a,b);
    double t=t1+t2;
    return t;
  }

  
  
  void printOrbitals(int num=-1) const
  {
    if(num==-1)
      num=HFOrbitals.size();
    cout<<"num"<<setw(4)<<"\t"<<setw(4)<<"l"<<"\t"<<setw(4)<<"jj"<<"\t"<<setw(4)<<"tz"<<"\t"<<setw(10)<<"e"<<"\t"<<setw(5)<<"OccuP"<<endl;
    for(int i=0;i<num;i++)
      {
	double norm=0;
	for(int j=0;j<HFOrbitals[i].Coeff.size();j++)
	  {
	    norm+=HFOrbitals[i].Coeff[j]*HFOrbitals[i].Coeff[j];
	  }
      	cout<<HFOrbitals[i].l<<"\t"<<setw(4)<<HFOrbitals[i].jj<<"\t"<<setw(4)<<HFOrbitals[i].tz<<"\t"<<setw(10)<<HFOrbitals[i].e<<"\t"<<setw(5)<<HFOrbitals[i].OccuP<<"\t"<<norm<<endl;
      }
  }
  
  ///initialize
  void init();

  ///sort and fill orbitals.
  void fillOrbitals();

  ///construct matrix and diagonalize, generate new HF orbitals
  void diag();
  
  ///iteration
  void iter();

  ///cal. HF energy
  void calEnergy();

  vector<DataType> TwoBodyMat;
  
  int A,Z;
  vector< SHF_Orbital<DataType> > HFOrbitals;
  DataType Energy;
 private:
  int numIter;
  int count;
  const System_TzType * pSystem_Tz;
};
#endif
